library(terra)
setwd("C:/Users/Admin_Geoteca/Documents/Projet/CI/")Manipuler les rasters avec R et le package Terra
AfroMap’R 2025
Note
Ce support est inspiré par les informations de ces différentes sources : Géomatique avec R (Giraud et Pecout 2024), Analyse d’image raster (et télédétection) (Madelin 2021) et Image processing and all things raster (Nowosad 2021). Les données ont été adaptées au contexte ivoirien.
Documents indispensables : Cheat sheet r-base
Packages utilisés
terra: importer, manipuler et exporter des fichiers rasters.sf: importer, manipuler et exporter des données géographiques vectorielles.
Activer le package terra
Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.
install.packages("terra")Installer le package à partir de l’interface rstudio est bien sûr possible.
Importer des données avec terra
Importer des objets vecteurs
La ligne qui servira au transect est un objet vectoriel. On peut la charger en tant que tel avec la fonction vect() de terra. L’objet chargé est de classe SpatVector.
Il est aussi possible de la charger avec le package sf conçu pour manipuler les objets vectoriels en priorité. Pour cela, voir le module de l’école thématique portant sur les vecteurs : manipuler des objets vectoriels avec sf
Ici nous chargeons un objet vectoriel au format .shp
# Import du fichier de ligne
ligne <- vect("Data/shp/ligne.shp")
class(ligne)[1] "SpatVector"
attr(,"package")
[1] "terra"Il est aussi possible de charger des objets au format .gpkg. Nous chargeons le masque qui servira aux prochaines opérations.
# Charger un gpkg
#zone_etude <- vect("Data/shp/limites.gpkg")
zone_etude <- vect("Data/shp/lac_kossou2.gpkg")Importer des objets raster
La fonction pour importer des objets raster est rast() quelque soit le type d’image.
Importer des MNT
# MNT
MNT_nord <- rast("Data/DEM/DEM_N.tif")
MNT_sud <- rast("Data/DEM/DEM_S.tif")Importer des images satellites en bandes séparées
# Importer les différentes bandes des images sentinels
#La bande bleue
bleu <- rast("Data/S2/Sentinel2_07_01_2023_B02.tif")
#La bande verte
vert <- rast("Data/S2/Sentinel2_07_01_2023_B03.tif")
#La bande rouge
rouge <- rast("Data/S2/Sentinel2_07_01_2023_B04.tif")
#La bande proche infrarouge
proche_infrarouge <- rast("Data/S2/Sentinel2_07_01_2023_B08.tif")Les données sont chargées. À partir de ces bandes, nous pouvons créer une dalle ou une tuile dans le visible.
#On peut créer un objet qui contient déjà les trois bandes dans le visible
image_RGB = c(bleu, vert, rouge)
#Affichage des métadonnées
image_RGBclass : SpatRaster
dimensions : 20976, 10980, 3 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 799980, 909780, 690240, 9e+05 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 29N (EPSG:32629)
sources : Sentinel2_07_01_2023_B02.tif
Sentinel2_07_01_2023_B03.tif
Sentinel2_07_01_2023_B04.tif
names : T29NRJ_202~104329_B02, T29NRJ_202~104329_B03, T29NRJ_202~104329_B04
min values : 2321, 2141, 1935
max values : 4050, 4240, 5092 Détail des métadonnées
| Type d’information | Informations | Descriptions |
|---|---|---|
| class : | SpatRaster | Type d’objet |
| dimensions : | 20976, 10980, 3 (nrow, ncol, nlyr) | Le raster fait 20976 par 10980 pixels sur 3 couches |
| resolution : | 10, 10 (x, y) | Un pixel équivaut à 10 mètres par 10 mètres |
| extent : | 799980, 909780, 690240, 900000 (xmin, xmax, ymin, ymax) | Indique les coordonnées max et min en x et y |
| coord. ref. : | WGS 84 / UTM zone 29N (EPSG:32629) | Le système de coordonnées géographique sources |
| sources : | Sentinel2_07_01_2023_B02.tif | Nom des fichiers d’origine |
| Sentinel2_07_01_2023_B03.tif | ||
| Sentinel2_07_01_2023_B04.tif | ||
| names : | B02, B03, B04 | Le nom de chacune des bandes |
| min values : | 2321, 2141, 1935 | Les valeurs minimales par bandes |
| max values : | 4050, 4240, 5092 | Les valeurs maximales par bandes |
Importation d’un raster avec plusieurs bandes incluses
RGB_PIR <- rast("Data/S2/RGB_PIR.tif")
#Affichage des métadonnées
RGB_PIRclass : SpatRaster
dimensions : 20976, 10980, 4 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 799980, 909780, 690240, 9e+05 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 29N (EPSG:32629)
source : RGB_PIR.tif
names : B02, B03, B04, B08
min values : 2321, 2141, 1935, 1873
max values : 4050, 4240, 5092, 7456 Visualiser des images avec terra
Visualiser un raster
On peut faire varier l’affichage en paramétrant l’argument type = de la fonction plot(). Chaque type fait appel à d’autres arguments qui lui sont propres.
Contininous est le mode par défaut, range permet de définir une plage d’affichage, ici, entre les valeurs minimales et maximales.
Comme dans les logiciels de SIG, minmax() doit calculer les statistiques pour certaines taches, il suffit de rajouter le paramètre minmax(image, compute = TRUE),
# Avec le type continuous
plot(MNT_nord, type = "continuous", range = c(min(minmax(MNT_nord, compute = TRUE)), max(minmax(MNT_nord, compute = TRUE)))) 
Interval est pratique pour faire des classes rapides, le paramètre breaks est le nombre de classes voulues
# avec le type interval
plot(MNT_nord, type = "interval", breaks = 10)
classes permet d’afficher une classification. Si les données du raster sont déjà classées, alors il suffit d’indiquer type = 'classes'. Sinon il faut construire cette classification.
Dans notre exemple nous construisons la classification à la volée avec la fonction classify(). Ici, la classification se fait de dix en dix en commençant de zéro à la valeur maximale (pour visualiser la séquence n’hétitez pas à exécuter la commande seq(0, max(minmax(MNT_nord)), by = 10) dans la console).
# Avec le type = classes et la fonction classify
plot(x = classify(MNT_nord, seq(min(minmax(MNT_nord)), max(minmax(MNT_nord)), by = 10)), type = "classes")
|---------|---------|---------|---------|
=========================================
Visualiser des images satellites
Visualiser une bande seule
On peut visualiser les bandes une par une avec la fonction plot() en r-base.
Depuis lobjetbleu` contenant la bande bleue
# La bande unique contenant le bleu
plot(bleu)
Depuis l’objet RGB_PIR contenant toutes les bandes. Ici on affiche la bande verte en la sélectionnant avec son nom B03
# La bande nommée B03 dans l'image à quatre canaux
plot(RGB_PIR$B03)
Depuis l’objet RGB_PIR on affiche la bande proche infrarouge avec son index
# Le canal numéro 4 de l'image à 4 bande, soit le proche infrarouge ('col = grey(0:100/100)' permet de mettre la bande en noir et blanc)
plot(RGB_PIR[[4]], col = grey(0:100/100))
Visualiser plusieurs bandes
Pour visualiser plusieurs bandes en même temps, on utilise la fonction plotRGB() du package terra.
#On essaye d'afficher l'image
plotRGB(image_RGB)Error in grDevices::rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, maxColorValue = scale): l'intensité de couleur 1.01979, hors intervalle [0,1]La visualisation d’un raster multibande nécessite de connaître nos images.
Ce message d’erreur est dû au fait que les valeurs des pixels ne sont pas dans une plage de 0 à 255.
Connaître les valeurs min et max avec la fonction minmax() permet de renseigner la plage d’affichage des pixels.
valeurs_min_et_max <- minmax(image_RGB)
max(valeurs_min_et_max[2,])[1] 5092On peut renseigner cette valeur dans l’argument scale = de la fonction plotRGB(). L’ajout du paramètre stretch = "lin" (pour linear) permet d’étirer les valeurs (entre les extrêmes) et d’améliorer le contraste de l’image.
plotRGB(image_RGB, scale = max(valeurs_min_et_max[2,]), stretch="lin")
Il peut être nécessaire de préciser à quelles couleurs correspondent les bandes. Ici il s’agit de l’index des couches dans l’objet ìmage_RGB, mais on peut utiliser leur nom.
#Affichage à partir de l'intex
plotRGB(image_RGB, r = 3, g = 2, b = 1, scale = max(valeurs_min_et_max[2,]), stretch="lin")
Point déchargement de la ram pour soulager nos ordinateurs
Les images peuvent être lourdes pour nos ordinateurs. Il est possible de supprimer les objets de la ram avec la fonction rm().
# Nous n'aurons plus besoin de cet objet à l'avenir
rm(image_RGB)Visualiser une image multibande depuis des objets distincts
# Pour afficher une image en créant un raster multibande pour la commande même
plotRGB(x = c(rouge, vert, bleu), scale = max(minmax(c(bleu, vert, rouge))[2,]), stretch = "lin")
Visualiser une image multibande depuis un objet contenant toutes les bandes.
Ici on précise la couleur de chaque bande pour obtenir une image en vraies couleurs.
# Pour afficher le spectre du visible à partir de l'image multicanal. On définit les couleurs que l'on veut attribuer à chaque nom de bande.
plotRGB(RGB_PIR, r = "B04", g="B03", b="B02", scale= max(minmax(RGB_PIR)[2,]), stretch="lin")
Ici en fausse couleur, qui permet de mettre en valeur l’activité chlorophyllienne, la végétation (pixels vers le rouge) :
# En fausse couleur
plotRGB(RGB_PIR, r = "B08", g="B04", b="B03", scale = max(minmax(RGB_PIR)[2,]), stretch="lin")
Point déchargement de la ram pour soulager nos ordinateurs
Les images peuvent être lourdes pour nos ordinateurs. Il est possible de supprimer les objets de la ram avec la fonction rm().
# Nous n'aurons plus besoin de ces objets à l'avenir
#rm(bleu, rouge)Fusionner, reprojeter et découper des dalles
L’objectif de cet exercice est de parvenir à la réalisation d’un profil topographique qui traverse le lac de Kossou à partir des deux tuiles MNT.
Fusionner deux tuiles
Nous allons fusionner les deux tuiles chargées dans la Section 1.2.1
par(mfrow = c(2, 1)) # permet de jouxter deux figures (une ligne, deux colonnes).
plot(MNT_nord)
plot(MNT_sud)
Dans un script r il faut utiliser la fonction dev.off() pour remettre les paramètres d’affichage à zéro. Ici ceci permet de revenir à l’affichage d’une seule figure à la fois.
Pour fusionner deux tuiles, on utilise la fonction merge()
# La fusion est possible avec la fonction merge.
MNT <- merge(MNT_nord, MNT_sud)
|---------|---------|---------|---------|
=========================================
# Affichage des métadonnées
MNTclass : SpatRaster
dimensions : 23633, 22741, 1 (nrow, ncol, nlyr)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : -8.745139, -2.428194, 4.185694, 10.75042 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : spat_a6486aab7431_42568_zxzxx9jyyx6UnnI.tif
varname : DEM_N
name : DEM_N
min value : -22.97385
max value : 1755.34058
Point déchargement de la ram pour soulager nos ordinateurs
Les images peuvent être lourdes pour nos ordinateurs. Il est possible de supprimer les objets de la ram avec la fonction rm().
# Nous n'aurons plus besoin de ces objets à l'avenir
rm(MNT_nord, MNT_sud)Regardons le résultat
# Regardons le résultat
plot(MNT)
Agréger, ou comment baisser la résolution d’un raster
Les images peuvent être lourdes pour nos machines, pour de nombreuses raisons, il peut être utile de baisser la résolution de l’image. Pour cela, la fonction aggregate() va nous aider. Le paramètre fact = n renseigne le nombre de pixels que l’on veut en horizontal/vertival. fun = permet de préciser la logique par lequel vont être agrégées les valeurs des pixels (mean, max, min, median, sum, etc.)
#Changeons la résolution du MNT avec un groupement de 100 pixels par 100 pixels
MNT_allege <- aggregate(MNT, fact = 100, fun = 'mean')
plot(MNT_allege)
Désagréger ou augmenter la résolution d’un raster
Inversement, lorsque des traitements impliquent plusieurs rasters, il peut être nécessaire que les rasters soient à la même résolution. Il est possible d’augmenter la résolution avec la fonction disagg(). Comme la fonction aggregate(), il faut renseigner le facteur, le nombre de pixels voulus par pixel en horizontal/vertical. method renseigne la méthode de désagrégation (near = plus proche voisin ou bilinear).
#Changeons la résolution pour 2 pixels par 2 pixels pour chaque pixel
MNT_enrichi <- disagg(MNT, fact = 2, method = 'bilinear')Pour connaître la résolution des pixels, il est possible de regarder les métadonnées ou d’utiliser directement la fonction res(). L’unité est celle utilisée dans le système de coordonnées.
#Résolution des MNT
res(MNT)[1] 0.0002777778 0.0002777778res(MNT_allege)[1] 0.02777778 0.02777778#res(MNT_enrichi):::.callout-tip title=“Point déchargement de la ram pour soulager nos ordinateurs” collapse=“false”} Les images peuvent être lourdes pour nos ordinateurs. Il est possible de supprimer les objets de la ram avec la fonction rm().
# Nous n'aurons plus besoin de ces objets à l'avenir
rm(MNT_allege, MNT_enrichi):::
Projeter
Visualisons les deux objets qui serviront au transect, le MNT et la ligne vectorielle chargée plus haut dans la Section 1.1 .
#Affichons le MNT avec la ligne cette ligne
plot(MNT)
plot(ligne, add = T) # add = T permet d'ajouter la ligne à un plot déjà ouvert
La ligne ne s’affiche pas : vérifions si les projections sont identiques avec la fonction same.crs()
same.crs(MNT, ligne)[1] FALSELes objets n’ont pas la même projection, pour les connaître regardons leurs métadonnées
MNTclass : SpatRaster
dimensions : 23633, 22741, 1 (nrow, ncol, nlyr)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : -8.745139, -2.428194, 4.185694, 10.75042 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : spat_a6486aab7431_42568_zxzxx9jyyx6UnnI.tif
varname : DEM_N
name : DEM_N
min value : -22.97385
max value : 1755.34058 ligne class : SpatVector
geometry : lines
dimensions : 1, 1 (geometries, attributes)
extent : 858699.9, 909170.1, 772839.1, 825722.2 (xmin, xmax, ymin, ymax)
source : ligne.shp
coord. ref. : WGS 84 / UTM zone 29N (EPSG:32629)
names : id
type : <num>
values : NaNLe raster est en WGS 84 EPSG:4326, et la ligne en WGS 84 EPSG:32632.
Projetons le raster pour qu’il corresponde à la ligne avec la fonction project() et la fonction crs() permettant de récupérer le crs d’un objet.
MNT <- project(x = MNT, crs(ligne))
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Réessayons d'afficher la ligne sur le MNT
plot(MNT)
plot(ligne, add = T)
Découper un raster
Il est possible de découper un raster de deux façons, soit avec la fonction crop() qui prend en compte l’emprise (la bbox) d’un objet, soit avec la fonction mask() qui découpe précisément selon la forme de l’objet.
Pour notre exemple nous allons utiliser la fonction crop() pour découper le MNT à l’emprise de la ligne.
MNT_transect <- crop(MNT, ligne)
plot(MNT_transect)
plot(ligne, add = T)
On peut aussi découper un raster selon un masque chargé dans la Section 1.1. Pour cela on utilise la fonction mask()
MNT_mask <- mask(MNT, zone_etude)
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
plot(MNT_mask)
plot(zone_etude, add = TRUE, border = "red", lwd = 4) # 'border =' permet de colorer uniquement la bordure du polygone, lwd = correspond à l'épaisseur de ce contour.
Le zone_etude ne réduit pas l’emprise de la zone d’étude, mais remplace les valeurs en dehors du masque pas NA.
Pour obtenir une zone d’étude à la bonne emprise et découper, on combine les deux fonctions crop() et mask().
MNT_mask <- crop(MNT_mask, zone_etude)
plot(MNT_mask)
plot(zone_etude, add = TRUE, border = "red", lwd = 4)
Intersecter un vecteur et un raster
On peut extraire les valeurs d’un raster en fonction de la forme d’un vecteur.
Ici on extrait la valeur des cellules du raster se trouvant sur la ligne.
transect_points <- extract(x = MNT_transect, y = ligne)Pour savoir le type d’objet que nous avons, ‘class()’ permet de nous aider.
# Interrogation du type d'objet
class(transect_points)[1] "data.frame"L’objet récupéré est de type data.frame. La fonction head() permet d’afficher les premières lignes d’une data.frame et la fonction rnow() pour avoir le nombre de lignes.
# Visualisation des premières lignes de la df
head(transect_points) ID DEM_N
1 1 229.9107
2 1 230.0214
3 1 227.9192
4 1 228.5510
5 1 227.1072
6 1 224.9774#Nombre de ligne
nrow(transect_points)[1] 3362
Note utile : mesurer la longueur des différents objets
Suivant le type d’objet (connaissable avec la fonction classe()),
voici comment les mesurer :
ncol() : permet de connaître le nombre de colonnes dans une table
# Nombre de colonnes
ncol(MNT)[1] 22870ncol(transect_points)[1] 2nrow() : permet de connaître le nombre de lignes dans une table
#Nombre de ligne
nrow(MNT)[1] 23859nrow(transect_points)[1] 3362lenght() : permet de connaître le nombre d’objets dans un vecteur (c())
# Nombre d'objets dans un vecteur
length(transect_points$DEM_N)[1] 3362length(c(4, 8, 9, 4, 2, 3))[1] 6nlyr() : permet de connaître le nombre de couches que contient un SpatRaster
# Nombre de couches dans un objet SpatRaster
nlyr(RGB_PIR)[1] 4ncell() : permet de connaître le nombre de pixels dans un raster
# Nombre de pixels
ncell(MNT)[1] 545655330On a intersecté 3362 cellules.
Visualisons notre coupe topographique :
# Affichons le MNT et le profil sur un même plot
par(mfrow = c(2, 1)) # permet de jouxter deux figures (deux lignes, une colonne).
plot(MNT_transect)
plot(ligne, add = T)
title("Zone de la coupe")
plot(x = 1:nrow(transect_points),
y = transect_points$DEM_N,
type = 'line',
xlab = "Distance au début du transect",
ylab = "Altitude (en mètre)")
title("Profil topographique")
Réaliser la coupe avec SF
Réaliser la coupe avec sf permet de définir un nombre de points à intersecter.
# La méthode pour le profil consiste à créer des points réguliers le long de la ligne, d'extraire les valeurs de pixel du raster pour chaque point et créer un graphique. Nous allons utiliser le package sf pour créer les points et donc, de convertir l'objet SpatVector en objet sf.
library(sf)
# st_as_sf() pour la conversion vers un objet sf
transect_points <- st_as_sf(st_sample(x = st_as_sf(transect), type = "regular", size = 1000), crs = 32632)
# La fonction extract permet d'extraire les valeurs du raster
transect_points <- extract(x = MNT_transect, y = transect_points)
# Voyons à quoi ressemblent les données avec les premières lignes avec la fonction head()
head(transect_points)
# Affichons le profil avec plot
plot(x = transect_points$DEM_N, y = 1:1000)
# Il est nécessaire de préciser les deux axes et il est nécessaire que les deux séries aient le même nombre de valeurs. La fonction nrow() permet de connaître le nombre de lignes dans un tableau. Changeons les points par des lignes
plot(x = 1:nrow(transect_points), y = transect_points$DEM_N, type = 'line')
# Affichons le MNT et le profil sur un même plot
par(mfrow = c(2, 1)) # permet de jouxter deux figures (deux lignes, une colonne).
plot(MNT_transect)
plot(transect, add = T)
plot(x = 1:nrow(transect_points), y = transect_points$DEM_N, type = 'line')Algèbre spatiale
On peut manipuler les valeurs des rasters selon quatre grands types d’opérations d’algèbre spatiale :
Locale : chaque pixel indépendamment, par exemple calculer un NDVI
Focale : la valeur du pixel dépend de la valeur des pixels l’entourant
Globale et Zonale : la valeur du pixel dépend de l’ensemble ou d’une zone de son espace de référence. Lorsque l’opération est globale, il s’agit de l’emprise de l’image, lorsque l’opération est zonale, il s’agit d’une partie de l’emprise de l’image.
Opération globale
Les opérations globales vont interroger l’ensemble des valeurs d’un raster pour obtenir un résultat, il s’agit souvent d’opération statistique.
terra regroupe quelques opérations dans la fonction global()
Essayons de mieux connaître nos images
Essayons de connaître l’altitude moyenne, minimale, maximale
global(MNT, fun = 'mean', na.rm = TRUE) mean
DEM_N 244.749global(MNT, fun = 'max', na.rm = TRUE) max
DEM_N 1752.367global(MNT, fun = 'min', na.rm = TRUE) min
DEM_N -20.36014Essayons d’afficher un histogramme de notre MNT
hist(MNT)
En extra : Pourquoi pas essayer de créer des histogrammes manuellement?
Essayons d’afficher un histogramme de notre image satellite avec les quatre bandes
plot(y = c(freq(round(bleu, -2))$count), x = c(freq(round(bleu, -2))$value), type = "line", col = "blue", xlim = c(0,6000), ylab = "nombre de pixel", xlab = "Valeur des pixels")
lines(y = c(freq(round(vert, -2))$count), x = c(freq(round(vert, -2))$value), col = "green")
lines(y = c(freq(round(rouge, -2))$count), x = c(freq(round(rouge, -2))$value), col = "red")
lines(y = c(freq(round(proche_infrarouge, -2))$count), x = c(freq(round(proche_infrarouge, -2))$value), col = "purple") Opérations locales
Isolons les surfaces en eau du Sebkha Sidi El Hani grâce au NDWI.
Voici notre zone d’étude
plot(MNT_mask)
La première étape est de réduire l’emprise des bandes pour qu’elles correspondent à l’emprise de la zone d’étude.
Pour cela nous utilisons le crop :
# Nous aurons besoin de la bande 3 (vert), de la bande 8 (Proche infrarouge) et d'un masque pour alléger le calcul
kossou_vert <- crop(vert, zone_etude)
kossou_Proche_infrarouge <- crop(proche_infrarouge, zone_etude)
par(mfrow = c(1,2))
plot(kossou_vert, main = "Vert")
plot(kossou_Proche_infrarouge, main = "Proche infrarouge")
dev.off()null device
1 Puis on calcule le NDWI
L’indice d’eau par différence normalisée NDWI (Normalized Difference Water Index)
Le NDWI est un indice normalisé permettant de mettre en évidence les surfaces en eau sur une image satellite, les zones ± humides.
NDWI = (vert – proche infrarouge)/(vert + proche infrarouge)
Pour Sentinel : NDWI = (Band 3 – Band 8)/(Band 3 + Band 8)
# Calculons le NDVI sous le même principe de la calculatrice raster d'Arcmap ou Qgis
NDWI <- (kossou_vert - kossou_Proche_infrarouge)/(kossou_vert + kossou_Proche_infrarouge)
plot(NDWI)
Allons plus loin dans l’opération locale, créons un raster binaire entre surfaces humides/en eau et les autres
On crée un histogramme pour visualiser la répartition du NDWI, éventuellement faire un ou des seuillages.
# La fonction hist permet de créer un histogramme, ce qui nous permet de définir un seuil
hist(NDWI)
abline(v = 0)
On peut choisir de visualiser le résultat du NDWI de façon binaire : eau / non-eau
# Création d'un raster True / False : toutes les valeurs supérieures à zéro sont TRUE les autres FALSE
NDWI_binaire <- NDWI > 0
head(NDWI_binaire) T29NRJ_20230107T104329_B03
1 FALSE
2 FALSE
3 FALSE
4 FALSE
5 FALSE
6 FALSEplot(NDWI_binaire )
Mais on peut aussi faire une classification sur la base des plages standards
| Borne inférieure | Borne supérieure | Interprétation |
|---|---|---|
| 0,2 | 1 | Surface de l’eau |
| 0 | 0,2 | Haute humidité |
| -0,3 | 0 | Sécheresse modérée, surfaces non aqueuses |
| -1 | -0,3 | Sécheresse, surfaces non aqueuses |
NDWI_class <- classify(NDWI,
c(minmax(NDWI)[1,],
-0.3,
0,
0.2,
minmax(NDWI)[2,]))
head(NDWI_class) T29NRJ_20230107T104329_B03
1 (-0.3 - 0]
2 (-0.3 - 0]
3 (-0.3 - 0]
4 (-0.3 - 0]
5 (-0.3 - 0]
6 (-0.3 - 0]plot(NDWI_class, type = "classes", levels = c("Sécheresse, surfaces non aqueuses",
"Sécheresse modérée, surfaces non aqueuses",
"Haute humidité",
"Surface de l’eau"))
Opération focale
Les opérations focales sont des traitements qui, pour chaque pixel, utilisent les valeurs des pixels l’entourant. La création d’une carte de pente, les calculs de TPI (terrain position index), et tous les outils relatifs à l’hydrologie sont des opérations focales.
Certaines opérations classiques sur des rasters d’altitude (i.e. des MNT modèles numériques de terrain) existent déjà dans le package terra et sont accessibles avec la fonction terrain(), c’est le cas du calcul de pente, du TPI ou encore de l’ombrage.
On indique l’objet sur lequel effectuer l’opération focale, le type d’opération, l’unité de mesure et surtout le nombre de voisins directs à comparer.
Quelques exemples de calculs possibles :
par(mfrow= c(2,2))
# Pente
pente <- terrain(x = MNT_mask, v = "slope", unit = "radians", neighbors = 8)
plot(pente)
title("pente")
# Aspect
aspect <- terrain(MNT_mask, "aspect", unit="radians", neighbors=8)
plot(aspect)
title("aspect")
# TPI
###plot(terrain(MNT_mask, "TPI", neighbors=8))
###title("TPI")
Avec la fonction focal() nous pouvons paramétrer précisément notre opération. L’argument w= spécifie le nombre de pixels distants du pixel central, l’argument fun= prend la valeur d’une fonction de calcul statistique.
Par exemple, il est possible d’adoucir le profil d’élévation calculé dans la Section 4.
D’abord on lisse les valeurs du MNT :
# moyenne focale
MNT_transect_focal_9 <- focal(x= MNT_transect, w = 9, fun = "mean")
MNT_transect_focal_19 <- focal(x= MNT_transect, w = 19, fun = "mean")Puis on intersecte le MNT avec la ligne vectorielle :
# La fonction extract permet d'extraire les valeurs du raster
transect_points_focal_9 <- extract(x = MNT_transect_focal_9, y = ligne)
transect_points_focal_19 <- extract(x = MNT_transect_focal_19, y = ligne)Puis on compare le résultat :
par(mfrow = c(3, 1)) # permet de jouxter deux figures (deux lignes, une colonne).
# On reprend le profil original
plot(x = 1:nrow(transect_points), y = transect_points$DEM_N, type = 'line', ylab = NA)
title("Coupe initiale")
#Pour le plot, supprimons les valeurs na, les opérations focales ne peuvent calculer les pixels en bordure d'image
plot(x = 1:nrow(na.omit(transect_points_focal_9)), y = na.omit(transect_points_focal_9$focal_mean), type = 'line', ylab = NA)
title("Coupe adoucie : focal = 9")
plot(x = 1:nrow(na.omit(transect_points_focal_19)), y = na.omit(transect_points_focal_19$focal_mean), type = 'line', ylab = NA, xlab = "Distance au point d'origine")
title("Coupe adoucie : focal = 19")
dev.off()null device
1 Représenter le terrain grâce aux opérations focales et locales
On peut combiner les objets construits grâce aux opérations focales et locales précédentes. Ici on construit une représentation du terrain mêlant l’ombrage et le NDWI.
Pour construire l’ombrage on mobilise nos objets pente et aspect dans la fonction shade() à laquelle on rajoute quelques paramètres.
ombrage <- shade(pente, aspect, angle = 40, direction = 300, normalize = FALSE)# ombrage
plot(ombrage, col = grey(0:100/100))
# NDWI
plot(NDWI_binaire, type = "classe", col = c(rgb(255, 255, 255, alpha = 0, maxColorValue = 255), "skyblue2"), add = T)
# MNT
plot(MNT_mask, alpha = 0.5, add = T)
En extra : Testons avec plusieurs images
Cet exemple permet de regarder l’évolution de la surface en eau à partir d’images prises à plusieurs dates. Bien que le code fasse appel à une boucle, il n’est pas bien compliqué puisqu’il reprend ce que l’on a vu jusqu’ici.
D’abord on charge les données automatiquement. Ici on fait la liste de tous les fichiers contenant les caractères “B03” dans le dossier Data/Sentinel2 avec la fonction list.files(), on charge les objets de cette liste avec la fonction rast() puis on les découpe avec crop() avec le masque fait plus haut.
On utilise la syntaxe avec les “pipe” |> pour indiquer une suite d’opération, mais on peut aussi les imbriquer.
On répète ensuite l’opération avec “B08”.
# Chargement et préparation de la bande du proche infrarouge (B08) suivant la logique des 'pipe'
#B08 <- list.files("Data/Sentinel2", pattern = "B08", full.names = T) |>
# rast() |>
# crop(y = zone_etude)
#
## Chargement et préparation de la bande du proche infrarouge (B08) suivant la logique des fonctions imbriquées
#B03 <- crop(rast(list.files("Data/Sentinel2", pattern = "B03", full.names = T)), zone_etude)Puis on réalise l’affichage
# Prépare une figure avec 5 lignes et 3 colonnes
#par(mfrow = c(3,3))
#
#for(i in 1:nlyr(B03)){ #1:nlyr(B03) permet de créer une suite de nombre commençant de 1 jusqu'au nombre de couches que contient l'objet B03. L'objet `i` aura pour #valeur chaque chiffre de la séquence 1:nlyr(B03)
# # Calcul le NDWI pour la couche i
# NDWI_i <- (B03[[i]] - B08[[i]]) / (B03[[i]] + B08[[i]])
# # Affiche le NDWI pour i
# plot(NDWI_i)
# # Dessine l'histogramme pour i
# hist(NDWI_i)
# # Affiche le raster binaire eau oui / non
# plot(classify(NDWI_i, c(minmax(NDWI_i)[1,], 0, minmax(NDWI_i)[2,])))
# # efface l'objet
# rm("NDWI_i")
#}Opérations zonales
Les opérations zonales sont similaires aux opérations globales. Cependant, au lieu de calculer des statistiques pour l’ensemble d’un raster, les calculs seront faits pour des emprises spatiales définies.
Avec un MNT et des limites administratives, il est possible de calculer pour chaque entité administrative, les altitudes moyennes, max, min, etc. Mais les possibilités sont plus nombreuses avec d’autres types de raster.
Essayons d’obtenir les altitudes moyennes pour chaque délégation de l’emprise de notre MNT.
Chargeons nos limites administratives
##Chargeons les limites avec la fonction vect()
#delegation <- vect("Data/geom/tun_admin.gpkg", layer = "delegation")
## Vérifion nos métadonnées
#delegationVérifions si le système de coordonnées est identique avec notre MNT avec same.crs()
# Les crs sont-ils identiques?
#same.crs(MNT, delegation) Reprojetons nos délégations
#delegation <- project(delegation, crs(MNT))
#plot(delegation)
#plot(MNT, add = T)Découpons nos limites administratives selon l’emprise du MNT
Créons un polygone qui épouse les bordures de notre MNT avec as.polygone()
MNT != -99999 permet de renseigner que nous voulons que les valeurs différentes de -99999, une altitude non réelle, ça fonctionne et c’est rapide.
#emprise <- as.polygons(MNT != -99999)
#delegation_decoupe <- crop(delegation, emprise)
#plot(MNT)
#plot(emprise, border = "red", lwd = 6, add = TRUE)
#plot(delegation_decoupe, add = T, lwd = 4)Créons un objet raster avec l’altitude moyenne avec la fonction zonale() et affichons le résultat
#zonalRast <- zonal(x = MNT, z = delegation_decoupe, fun = "mean", as.raster = TRUE)
#head(zonalRast)
#plot(zonalRast)Faisons la même chose, mais en créant des objets vecteurs
#zonalVect <- zonal(x = MNT, z = delegation_decoupe, fun = "mean", as.polygons = TRUE)
#head(zonalVect)
#summary(zonalVect)
#plot(zonalVect, col=rev(heat.colors(10)))Exporter un fichier
Vers .gpkg, vers .tif, vers format sf
r ne crée pas automatiquement les dossiers : il faut les créer en amont
Lorsqu’on enregistre un fichier, si le dossier n’existe pas, r ne va pas le créer par défaut. Exemple de code pour vérifier l’existence du dossier : dir.exists() informe par vrai ou faux si un répertoire existe. Si c’est vrai, il exécutera le code qui suit entre les accolades. Si on veut exécuter quelque chose si la condition est fausse, il suffit d’écrire else et écrire ce qu’il faut exécuter si faux.
#if(dir.exists("Data/Export/")){ # si le dosser Data/Export/ existe ...
# print("Le dossier existe déjà") # il faut afficher ce message : Le dossier existe déjà
# }else{ # sinon ...
# dir.create("Data/Export/") # on crée le dossier Data/Export/ ...
# print("Le dossier vient d'être créé") # Et on affiche le message Le dossier vient d'être créé
#}Exporter un raster avec writeRaster()
#writeRaster(MNT_transect_focal_19, "Data/Export/MNT_lisse.tif", overwrite = TRUE)Exporter un vecteur avec writeVector()
#writeVector(transect, "Data/Export/Ligne.shp", overwrite = TRUE)
#writeVector(transect, "Data/Export/Ligne.gpkg", overwrite = TRUE)Convertir un objet SpatVector vers un objet sf
#Vérifions que notre fichier est bien un objet SpatVector avec la fonction class()
#class(delegation_decoupe)Utilisons la fonction st_as_sf() pour convertir le type d’objet et vérifions avec la fonction class() :
#library(sf)
#delegation_decoupe_SF <- st_as_sf(delegation_decoupe)
#class(delegation_decoupe_SF)
#delegation_decoupe_SFRegardons si on trace notre MNT avec nos deux objets vecteurs, l’un sf, l’autre SpatVector, s’ils s’affichent et se superposent bien :
#Affichage du MNT en nuances de gris
#plot(MNT, col = grey(0:100/100))
##Affichage de notre objet sf en rouge
#plot(delegation_decoupe_SF, add = TRUE, border = "red", lwd = 10)
##Affichons notre SpatVector en blanc
#plot(delegation_decoupe, add = TRUE, lwd = 4, border = "white")
A propos de ce document
Ce support a été créé pour la semaine de formation franco-tunisienne GEO UNIV’R Tunisie 2024 - “Enseigner la statistique, la cartographie et l’analyse spatiale avec R qui se tient à Sousse en mai 2024.
Références
sessionInfo()R version 4.4.3 (2025-02-28 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.utf8 LC_CTYPE=French_France.utf8
[3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C
[5] LC_TIME=French_France.utf8
time zone: Europe/Paris
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] terra_1.8-21
loaded via a namespace (and not attached):
[1] digest_0.6.37 codetools_0.2-20 fastmap_1.2.0 xfun_0.50
[5] knitr_1.49 htmltools_0.5.8.1 rmarkdown_2.29 cli_3.6.3
[9] compiler_4.4.3 rstudioapi_0.17.1 tools_4.4.3 evaluate_1.0.3
[13] Rcpp_1.0.14 yaml_2.3.10 rlang_1.1.5 jsonlite_1.8.9
[17] htmlwidgets_1.6.4
Giraud, Timothée, et Hugues Pecout. 2024. Géomatique avec R. https://doi.org/https://doi.org/10.5281/zenodo.5906212.
Madelin, Malika. 2021. « Analyse d’images raster (et télédétection) ». https://mmadelin.github.io/sigr2021/SIGR2021_raster_MM.html.
Nowosad, Jakub. 2021. « Image processing and all things raster ». https://nowosad.github.io/SIGR2021/workshop2/workshop2.html.